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ABSTRACT 

In this paper we show how to use data on Lyman-limit and Damped Lyman-a 
absorption systems to derive the hydrogen ionization fractions and the distribu- 
tion of the face-on total gas column density. We consider axially symmetric, 
randomly oriented absorbers, ionized by an external background radiation field 
in order to relate the face-on total gas distribution to that of the neutral hydrogen 
observed along the line of sight. We devise a statistical procedure based on the 
Maximum Likelihood criterion, that is able to treat simultaneously data coming 
from different surveys and statistically recovers the "true" column densities in 
the presence of large uncertainties: this is especially important for Lyman-limit 
systems which leave an unmeasurable residual flux at wavelengths shorter than 
the Lyman break. We make use of simulated data to look for possible observa- 
tional biases and extensively test our procedure. For a large statistical sample of 
real data in the redshift range [1.75,3.25] (collected from all published surveys) 
our Maximum Likelihood procedure gives a power-law slope for the total hydro- 
gen distribution of —2.7. All together Lyman-limit systems therefore contain 
more gas than Damped Lyman-a systems. Analysis of data at other redshifts 
shows that more observations are needed to reach a compelling evidence for a 
cosmological evolution of the slope of the gas distribution. 

Subject headings: methods: statistical — quasars: absorption lines — catalogs 
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1. Introduction 

Absorption features in quasar spectra represent a powerful tool to investigate the for- 
mation and evolution of gaseous structures in the early Universe (Rauch 1998, and references 
therein), the efficiency of processes such as merging (Kauffmann 1996), gas depletion due to 
star formation (Wolfe et al. 1995; Storrie-Lombardi and Wolfe 2000), and ionization by the 
ultraviolet background radiation field (Weinberg et al. 1997; Savaglio et al. 1997). Also the 
level of cloud clustering may be investigated analyzing absorption spectra (Cristiani et al. 
1997). Lyman-limit absorption systems (hereafter LLS) are cosmological structures which 
contain enough neutral hydrogen to absorb Lyman continuum photons and produce a break 
in the QSO continuum flux level. This feature is detectable in moderate-resolution spectra 
whenever the opacity to Lyman continuum photons satisfies the condition tll ^ 1, equiva- 
lent to a neutral hydrogen column density iVjjj > 1.6 x 10 17 cm -2 . The column density of LLS 
can be determined whenever a residual flux beyond the Lyman break is detectable. Systems 
with iVjjj > 5x 10 19 cm -2 are also visible in moderate-resolution spectra since they give rise 
to Damped Lyman-a absorption lines with a reference frame equivalent width > 5 A. In 
this paper we shall refer to these as Damped Lyman-a absorption systems (hereafter DLS). 
The Hi column density of absorbers is usually very poorly determined observationally when 
5 x 10 17 < iVjjj < 5 x 10 19 cm~ 2 since the Lyman break is saturated, while the Lyman-a 
absorption line is not yet damped. Previous studies of absorbers with poorly determined 
column densities have been carried out using a coarse binning in iVjjj. But a coarse binning 
suffers by statistical problems, related to the fact that different choices of binning may lead to 
different results (e.g. Stengler-Larrea et al. (1995) for problems related to redshift binning), 
as well as by some negative consequences on the physical description, because it does not 
differentiate Hi column densities which correspond to different ionization conditions. This 
is particularly important for 10 17 < < 10 20 cm -2 since in this range a smoothly varying 
total gas column density distribution leads to a rapidly varying Hi column density distri- 
bution due to the rapid change of the H ionization fraction (Corbelli et al. 2000, hereafter 
Paper I). 

To tackle this problem, we have collected data on LLS and DLS from all the available 
literature and perform a comprehensive statistical analysis, based on a Maximum Likelihood 
procedure, which allows a joined fit of observations with different iVjjj sensitivities and 
uncertainties, with different ionization conditions and at different redshifts. 

The outline of this paper is the following: in Section 2 we discuss how the observed Hi 
column density distribution is related to the total hydrogen distribution for face-on absorbers. 
A set of analytical solutions for the cross section of non-spherical absorbers, randomly ori- 
ented in space, is reported in Appendix A. Section 3 contains a description of the database. 
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Section 4 details the statistical algorithm that we use to best fit the Hi column density distri- 
bution of absorbers. Section 4 also summarizes results relative to our database. In Section 5 
we perform realistic simulations of the data, which are used to test the effectiveness and the 
robustness of the statistical algorithm. 



2. Column densities relationships 

One physically meaningful quantity for investigating the formation and evolution of gas 
condensations is the distribution of their gaseous mass. However observations measure only 
the distribution of the neutral hydrogen column density along the line of sight, which may 
present a rather complex behavior even if the face-on total column density distribution is a 
power law. In this section we show how to transform one distribution into the other. The 
relationship between the total and the neutral column density has been investigated in detail 
in Paper I. Here we focus on the effect of randomly oriented absorbers which are either flat 
or of ellipsoidal shape. 



2.1. From the total to the neutral column density distribution 

We show first how a face-on total hydrogen column density distribution g±{Nn ± ) trans- 
lates into a face-on neutral hydrogen distribution /±(A^jj_|_). For a generic relationship 

-^ELl (-Whll) one can wr h e: 

h (JVhu) = 9± (JVki) (1) 

dl\ m± 

If g± is a pure power-law: 

9± (Wh±) = KN£ (2) 
f± (Whll) 1S a l so a power law when Nn ± = QN^ I± : 

f± (Nm±) = (3KQ-^N-f ± , (3) 

with £ = 1 + (3 (a — 1). For f3 < 1 (an ionization fraction monotonically decreasing with 
AfjjjJ and a > 1, f± (iVjjij_) is always shallower than g± (N^jJ. 

For a generic iVjjj_ (iVjjij.) relationship, /±(-/Vjj i± ) is not tied to a power law. For example 
when there is a sharp break in the slope of Nn ± ( A^jji_l ) a t Nb, f± is not continuous at JVj,, 
being larger on the side where its slope is steeper: 

f± (N m± ) = k 1 ( (a - 1) (N m± /N b y^e (N b - N m± ) + 
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(6 - 1) (N m± /N b y^Q (iV Hl± -N b )), (4) 

where 0(x) is the Heaviside step function. 

Similarly we can consider multiple breaks: for example if N^ l± = iVjj ± for large iVjj I± 
and iVjjjj_ = eiVji_|_ for small N-r i± (as for a constant ionization fraction), f± follows a power 
law with index —a in both regions, but the low-column density distribution is shifted by a 
factor e"" 1 with respect to the extrapolation of the high-column density one. 



2.2. Orientation-averaged differential cross section 

In this section we consider the effects of the random orientation of non-spherical ab- 
sorbers on the observed column density distribution. We derive a simple analytic solution 
for "infinitely thin homogeneous slabs" , namely those whose thickness (2h) is much smaller 
than their radial extension (2R), and whose face-on gas column density is constant. We then 
test a simple approximation for treating "finite thickness slabs" , using the exact solution for 
homogeneous ellipsoids, given in Appendix A. 

Let c(A / jj I ,//) be the differential cross section of an absorber, whose symmetry axis is 
tilted by an angle 9 with respect to the line of sight. Its differential cross section, averaged 
over orientations, is: 

1 f 1 

&(N m ) = - J a{N m ,n)dn (5) 

(where we have defined \x = cos 9). For infinitely thin slabs of surface area A and face-on 
column density Ajji±, the averaged cross section reads: 

1 P 1 AN 2 
*(Nm) = g J / H *(^Hi - AWH) dfi = -j^ 6 (N m - N m± ) , (6) 

where S(x) is the Dirac distribution. Therefore the averaged cross section vanishes for Ajjj < 
Ajjj_|_, and is proportional to Nj^ for Ajj : > Nji 1± . 

For absorbers of finite thickness the exact behavior of a depends on their detailed shape, 
but the following approximate formula holds: 

AN 2 

HNm) = -ttP © (Nm - JVhll) © (Nmj.R/h - N m ) . (7) 

iV Hl 

It is shown in Appendix A that the accuracy of Eq. (7) decreases as h/R grows larger. 
When h ~ R, a for ellipsoidal absorbers differs noticeably from Eq. (7); but, once Eq. (7) is 
convolved with realistic f± distributions, results keep close to the exact ones. 
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2.3. From the face-on to the line of sight A^jji distribution 

The line of sight iVjjj distribution, /(-/Vjjj), can be written as: 

/ (N m ) = J ^j^f±(N m ±) dN m± , (8) 
which, using Eq. (7), gives: 

1 f N m 

f (Nm) = TTir / N m J ± (N m± ) dN m± . (9) 

iV Hl JN m h/R 

This equation shows that /(-/Vjjj) is sensitive to the values of f± over the whole range 
[N Rl h/R, N m ] (i.e. [0, N m ] for infinitely thin slabs). 

If there is a single break in f±, as given by Eq. (4), the corresponding f(Nm) is: 

/ (N m ) = ( fr^ 1 ( min (iVHl ' iYfc)3 " a ~ min ( N n h / R > ^) 3_a ) + 

|^iV« 2 (max (N m , N h f^ - max (N m h/R, N b ) 3 ^jj . (10) 

For A^jjj < iVf, and for N b R/h < iVjjj, / follows power laws with indices — £i and —£2 
respectively. In the intermediate iVjjj range the behavior of / is more complex: if f± is 
steeper at TVjjj < N b , its slope approaches —3 for N b < TVjjj <C N b R/h while it bends to 
connect to the power-law regime at iVm > N b R/h; in the opposite case / suddenly rises for 
TVjjj > N b , and then smoothly connects with the high column density power- law branch. 

Fig. 1 shows qualitatively the features present in /(-/Vjjj) for a realistic -^h±(-^Hi±) 
relationship (see Paper I): Nmi oc N^ I± for N b i < Nm± < N b2 , Nn ± oc Nm± elsewhere. We 
choose a = 1.5, h/R = 0.1, iV fel = 10 17 cm" 2 , N b2 = 10 20 cm- 2 , and (3 = 0.001 to emphasize 
the various regimes described above. Note that f(Nm) at low column densities is lower 
by a factor (N bl / 'N^ 1 '^^ 1 ^ with respect to the extrapolation of the distribution at high 
column densities (dashed line). At intermediate column densities a dip in the distribution 
appears, approaching a —3 slope towards the low column density side. 

3. DLS&LLS surveys and our data sample 

In this section we show which information should be extracted from the large number 
of DLS&LLS surveys available in order to build up a homogeneous database for a global 
statistical approach. 
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lglO N HI ( Cm_2 ) 

Fig. 1. — For a given face-on iV Hl distribution (dotted line) we show the projected distribu- 
tion /(-/Vjjj) (solid line). The original face-on total gas distribution is also plotted (dashed 
line). Power-law indices in regions where /(-/Vjjj) approaches a power law are indicated. We 
have used a realistic relationship between N-r ± and iVjj I± given explicitly in the text. 
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In a LLS the value of A^jji can be determined from the ratio between /, the residual 
Lyman continuum flux on the blue side of the break, and I , the unabsorbed continuum flux: 

iVm = 1.6 x 10 17 r LL cm- 2 = 1.6 x 10 17 ln(/ //) cm" 2 . (11) 

If the residual flux / is too small to be measured, only a lower limit to iVjjj can be derived: this 
usually happens for tll ~ 3. Since there are measurements sensitive to Lyman continuum 
optical depths as small as 0.4, we include in our database all LLS surveys which are sensitive 
to JV H i > 10 16 - 81 cm" 2 . 

We include in the compilation all available searches for absorption lines with a rest 
frame equivalent width W > 5 A, for which Ajii can be estimated from: 

N m = 1.88 x 10 18 W 2 cm" 2 . (12) 

A complication comes from the fact that spurious lines with large equivalent width may 
result from a blending of weaker lines with metal line systems or Lyman-a forest lines. Only 
some surveys have sufficient spectral resolution to test, by fitting a Voigt profile, whether 
lines with W > 5 A are really damped. In our database we include both the column density 
value determined directly from Eq. (12) and that derived from the Voigt profile fit, whenever 
this is available. By comparing these two values we find that there is a bias in the iVjjj values 
estimated from the equivalent width, and derive the following statistical correction: 

Igio N m , corr = 6.261 + 0.705 lg 10 N m , w . (13) 

We apply this correction to all absorption lines for which only the W value is available. If 
only a lower limit to Tll can be derived from the Lyman continuum absorption, searches 
for the corresponding damped Lyman-a line can be used to establish a value or an upper 
limit to the Hi column density. If the corresponding Lyman-a absorption line has not been 
searched for, we take an upper limit of 5 x 10 21 cm -2 , which is the highest measured value 
of A"jji in our database. 

Unfortunately upper limits to saturated LLS leave large uncertainties in the column 
density range 5 x 10 17 ^ A r jji ~ 5 x 10 19 cm~ 2 , where theoretical predictions are very sensitive 
to the detailed shape of the -/Vhj-C-^Hll) relationship. Hence it is essential to find a correct 
statistical approach for evaluating the iVjji distribution in the presence of large uncertainties 
to determine the ionization conditions of the gas as well as the total gas column density 
distribution in LLS and DLS. 

We have collected data from the following references: Tytler (1982); Bechtold et al. 
(1984); Wolfe et al. (1986); Tytler (1987); Lanzetta (1988); Sargent, Steidel and Boksenberg 
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(1989); Turnshek et al. (1989); Lanzetta (1991); Lanzetta et al. (1991); Courvoisier and 
Paltani (1992); Bahcall et al. (1993); Storrie-Lombardi et al. (1994); Lanzetta et al. (1995); 
Stengler-Larrea et al. (1995); Wolfe et al. (1995); Storrie-Lombardi et al. (1996); Jannuzi 
et al. (1998); Storrie-Lombardi and Wolfe (2000). In order to ensure an accurate estimate 
of the coverage of the survey we have excluded from the compilation the full spectrum or 
part of it whenever the continuum was visually found to be very noisy or non-uniform (e.g. 
some spectrum edges where the flux is below the 2-a level in Lanzetta et al. (1995) or the 
objects 1130-106Y, PKS1206+459, MC1215+113 in Bahcall et al. (1993)). The final number 
of directions i.e. of QSOs included is 661. We have decided to split each direction in parts 
(hereafter called "paths" ) such that each path is defined by homogeneous search parameters 
and contains at most one detection. In other words one direction may correspond to more 
than one path if different sensitivities are used, or if more than one absorber is present along 
the line of sight. Our database is available in electronic form upon request to the authors. 
Parameters contained in the data table can be divided into two groups: those describing the 
search, and those related to the detection. 

The main parameters associated with the search coverage are: (i) the redshift path 
limits, i.e. the upper and lower redshift for one direction in space observed with a given 
sensitivity; (ii) the sensitivity, specified by the threshold above which iVjjj was detectable. 

For each direction the upper redshift was set to 5000 km s -1 less than the QSO redshift, 
whenever the spectral coverage extends beyond it. We have usually taken as the lowest 
redshift that given by the observers. This does not coincide with the spectral coverage 
boundary in the presence of "shadowing effect" . The "shadowing effect" intervenes whenever 
the continuum flux, absorbed by a LLS at z*, gets below detectability: in this case the lower 
limit of the path has been set to z*. If no DLS searches have been performed no further paths 
are appended along that direction; otherwise, due to the different frequency of the Lyman-a 
line respect to the H ionization threshold, further paths are considered between z* and the 
Lyman-a shadowing redshift (z = 0.752* — 0.25). The lowest redshift of the LLS search in 
the absence of shadowing is reported in the database as well, since in Section 5 it will be 
used for simulating data samples. The sensitivity for a Lyman-limit search depends on the 
minimum opacity detectable in the survey: typically r LL = 0.4, 1.0 or 1.5, corresponding 
to column density thresholds \g 10 N th =16.81, 17.20 and 17.38 respectively. We have set to 
0.4 the minimum detectable tll for observations listed in Bahcall et al. (1993) and Tytler 
(1987); we have instead set to 1.0 the threshold for all the observations listed in Lanzetta 
(1991) and Sargent, Steidel and Boksenberg (1989) which were revised by Stengler-Larrea 
et al. (1995). For absorption line searches with rest frame equivalent width W > 5 A the 
sensitivity of the path was set to lg 10 = 20.13 (since we have applied the correction given 
in Eq. (13)). If both LLS and DLS searches were performed on the same path the threshold 
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of the path was set to the lowest of the two, except in the presence of shadowing effects. 

The parameters associated with the detections are: (i) the redshift of the detected 
absorber; (ii) the decimal logarithm of the estimated Hi column density of the absorber 
along the line of sight (lgio^m); i™) the minimum and the maximum allowed values for 
lgio-^Hi! fa) the Voigt correction factors. 

When lg 10 A^jji could not be determined from t ll or W measures, but only limits for it 
are available, we have reported as lg 10 iVjjj the arithmetic mean between the minimum and 
maximum allowed values. When instead the value of lg 10 Am could be directly determined we 
have estimated the minimum and maximum allowed values by assuming a standard spread. 
In the case of non-saturated LLS we have used a spread of ±0.2 in decimal logarithm, except 
for the cases listed in Jannuzi et al. (1998) for which we took the uncertainties quoted by 
the authors. In the case of DLS we have used again a spread of ±0.2 when the line profile 
has been resolved and the Voigt profile has been fitted; otherwise we have used ±0.4. The 
uncertainties in lg 10 iVjjj that we quote are usually larger than those that can be directly 
derived from uncertainties on either tll or W. This is because there are other sources of 
errors, such as the determination of the level of the continuum (which affects both W and r LL 
estimates), or the possible blending with other small absorption lines. The Voigt correction 
factors are the differences between the lg 10 iVjjj values derived from the Voigt-profile fit and 
those from W (Eq. (12)). 

Two further columns contain flags, which are different from whenever LLS and DLS 
searches respectively have been performed on a given path. A value of —1 implies that no 
systems were found, a value of 1 means that the column density estimate is quite accurate, 
2 that it is less accurate and 3 that it is very uncertain, as for DLS declared "non damped" . 
In this last case the column density value for the DLS has been set to zero, or used as an 
upper limit if saturated LLS are detected at the same redshift. 

The coverage of the whole sample, with and without shadowing, for various search 
sensitivities is shown in Fig. 2, where the inhomogeneity of the actual coverage is evident. 

4. Statistical Analysis of the Data 

Various authors have analyzed Lyman-a absorbers data to determine the Ajjj distri- 
bution and/or its evolution with time. Nevertheless no statistical analysis of the iVjjj dis- 
tribution in LLS and DLS has been attempted so far using all available data and without 
a pre-defined coarse binning. In the following we show how to perform a global unbinned 
statistical analysis of heterogeneous data. 
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Fig. 2. — The redshift coverage of our sample for the various thresholds in column density 
(lg 10 N th =16.81, 17.20, 17.38 and 20.13). The continuum lines refer to the actual (shadowed) 
coverage, while the dotted ones refer to the case without shadowing. 
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4.1. Maximum Likelihood Analysis 

Our approach for estimating the iVjjj distribution is similar to the Maximum Likelihood 
analysis used by Storrie-Lombardi et al. (1996), with the following two improvements: (i) we 
use a generic (not just a pure power-law) distribution for iVjjj; (it) we take into account 
uncertainties in the iVjjj measurements. In order to simplify the notation we shall use the 
symbol N instead of iVjji in this section as well as in Section 4.2. 

Let f(N, z) be the distribution of column density iV at redshift z. The average number 
of detections along a given direction is f(N, z) Sv, where Sv = 5N5z is the area of an in- 
finitesimal cell in the N-z plane. For infinitesimal cells the probability of having more than 
one detection in the same cell is negligible. The probabilities of having zero or one detection 
in a given cell are respectively: 

P (N, z) = exp (-/ Sv) ; P,(N, z) = (fSv) exp (-/ Sv) . (14) 



The likelihood for a given data set is the product of the probabilities over "all" cells 
in the N-z plane and over all directions in space. This product can be separated into two 
factors, one relative to the "empty" cells, and the other to the "full" cells: 

C =X\P= \{ P Qti \{ P hj = H exp(/,H J] (/». (15) 
i=all i=empty j=full i=all j=full 

We know the total number of full cells (i.e. the number of detections, hereafter p) but we 
might not be able to position the detections into specific cells in the N-z plane due to the 
large uncertainties of N. For instance, if a detection has Nj in the range [N min j, N maXt j], the 
absorber is in one of the cells of this N range, but we don't know in which one: the probability 
for this event is the sum of probabilities relative to the cells (the same considerations would 
apply to errors in z if redshifts are not well determined). In the limit of small cell sizes, 
the sum is then approximated by an integral of / over the allowed range of N, and the 
logarithmic value of C Q reads: 

w(N,z)f(N,z)dNdz+J2^\ f{N, Zj )dNSz\, (16) 

j = l \" Nminj J 

where w(N, z) is the number of different paths at a given z with a column density threshold 
smaller than N. The first term in Eq. (16) is related only to the parameters of the search, 
and can be evaluated using the normalization condition for /: 

J w(N, z)f(N, z) dNdz = p. (17) 
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We can rewrite Eq. (16) using the above normalization condition: 

p ( fN maxj \ 

\nC = p(\n(5z)-l)+J2 ln [ / f(N, Zj )dN\. (18) 

j = l y' N min,j J 

If column densities have been measured with a rather small uncertainty for a number 
Pi of absorbers, Eq. (18) can be approximated as: 

ln£ = p(ln(^)-l)+^ln(5iV J )+^ln(/(iV„ % ))+^ln / ' f(N, Zj ) dN , (19) 

J = l J = l J=P1 + 1 Y N min,j J 

which allows a faster numerical computation. The terms with SNj and Sz are constant and 
can be neglected while searching for the maximum value of C Q . The limiting case p± = p 
reduces Eq. (19) to the formula in Storrie-Lombardi et al. (1996). 



4.2. The fit goodness and the true column density distribution 

The Maximum Likelihood analysis does not contain any statistical test on the goodness 
of the "best" solution for /, since the absolute value of C Q is undetermined by a constant 
factor. Such a test can be done for example by comparing the theoretical versus the ob- 
servational cumulative function. For large uncertainties of the column density values, this 
comparison requires a reasonable guess at the "true" value of ./V for each individual detec- 
tion. A "reasonable" guess does not necessarily mean that each N value should reproduce 
the true one but that their overall distribution cannot be distinguished from the original 
one. We will show how to find this "reasonable" guess in the presence of large measurement 
uncertainties in N in order to avoid the use of a coarse binning. 

For a given absorber j the probability of having a true column density smaller than N 
in the allowed range [N min j, N maxJ } is: 

f(N', Zj )dN'/ f(N', Zj )dN'. (20) 

- ^ min,j J N min.j 

Rj behave as random quantities and should have a uniform distribution in the interval [0, 1], 
independently of the shape of the distribution function /. 

The probability of having a detection with column density smaller than iV in the overall 
iV range, is: 

/N . poo 

w(N', Zj )f(N', Zj )dN' / w(N', Zj )f(N', Zj )dN'. (21) 
-hh ' JN th 
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The distribution of Cj should also be uniform in the interval [0, 1], and is equivalent to the 
N cumulative distribution for a narrow range of redshift. 

For a sample of Xj (sorted for increasing Xj values, with < Xj < 1 and 1 < j < p) the 
cumulative distribution can be written as: 



S(x)= { j/p, Xj < x < x j+1 (22) 



o, 


X < X\ 


j/p, 


Xj ^ X <C Xj^-\ 


1, 


Xp < ~^~ X • 


given 


a distribution 



one is to use the Kolmogorov-Smirnov test on the maximum deviation between the two 
cumulative functions: 

D = max \S(x) - x\ ; (23) 

0<x<l 

the significance of D being given by the function Qks(\/pD) (Kendall and Stuart 1967). 

However if the maximum discrepancy between the two distributions originates from some 
bias, the Kolmogorov-Smirnov test becomes completely insensitive to the overall match. Such 
effect is reduced if instead of D we use the quantity U : 

U 2 =pJ^ (s(x)-x- (S(x') - x') dx'^j dx. (24) 

U varies with Xj more smoothly than D does, and its significance is Qks(^pU) (Kendall and 
Stuart 1967). We will use U for testing the uniformity of Rj and Cj since it is practically 
more effective and robust than D to ensure a good matching over the whole iV range. 

For any given distribution / we compute a "best guess" for the "true" column densities 
by maximizing the product of the significances relative to Ur and to Uc (hereafter referred 
to as Cr and Cc, respectively). There are always infinite guesses which allow a uniform 
distribution for Rj. However Rj and Cj are mutually dependent and their relationship 
depends on / (see Eqs. (20) and (21)). For an incorrect choice of / or of the estimated 
ranges [N min j, N max j] it may then become impossible to derive uniform distributions for Rj 
and Cj at the same time. 



4.3. Numerical results 

Given an original distribution function for the face-on total column density of absorbers 
of the form: 

9±(N RL ,z) = K(l + z)''N£ (25) 
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and a functional relationship between Nn ± and iVjj I± , a routine derives the Maximum Like- 
lihood values of a and 7 for an average absorber cross section a (Eq. (7)). K is fixed by 
the normalization condition to the total number of detections (Eq. (17)); 7 accounts both 
for cosmological effects, and for the physical evolution of the absorbers. Here we assume no 
evolution in the power-law index a in a given z interval and a pure power-law dependence 
for g±, but it would be straightforward to implement our routine if necessary to change these 
assumptions. Physically meaningful Ajj i± -7Vji ± profiles in various redshift ranges have been 
computed by simulating the ionization structure of gaseous slabs (see Paper I). To identify 
each profile we use the parameter X defined as: 

X = lg 10 (AW^Hi±) for N m = 1.6 x 10 17 cm- 2 . (26) 

(the corresponding iVjj I± depends on a) . An increasing value of X corresponds to a neutral- 
to- ionized transition getting sharper and occurring at higher iVjj. After determining the best 
a for fixed values of the parameters X, 7 and (h/R) we can determine the best values for 
these 3 parameters by comparing the relative Maximum Likelihood solutions. We use the 
Likelihood given in Eq. (19) with pi as the number of detections with an A^jji uncertainty, 
Igio N max — lg 10 N min < A (lg 10 iVjjj). Values for Rj and Cj are computed for all the remaining 
ip— Pi) detections. Such a computation is not relevant for the Likelihood maximization itself, 
but it is useful in order to compare the theoretical and the observational cumulative. Since 
it is not always possible to obtain uniformly distributed Rj and Cj for any set of model 
parameters, the uniformity level of Rj and Cj distributions is a further test on the quality 
of the fit to /(Vjji) distribution. For this reason the routine maximizes the product of Cr 
and Cc- 

The main results of our statistical fit procedure to data in the redshift range [1.75,3.25], 
where most of the observations are, have been discussed in Paper I, together with the relative 
cosmological implications. The Maximum Likelihood solution, a = 2.70 and X = 2.82, is 
obtained for 7 = 1.0 and h/R = 0.2 but a and X are only weakly sensitive to reasonable 7 
and h/R variations. We now investigate the possible dependence of a and of the cumulative 
function on the limiting value of A(lg 10 iVjj I ) for the treatment of errors, keeping X fixed to 
2.82. 

We consider 4 different cases: (a) no error treatment, i.e. for all absorbers we use the 
estimated Hi column density (either the measured value or the middle point in the allowed 
lg 10 iVjjj range); (b) and (c) error treatment has been applied to all the detections with 
A(lg 10 A^ Hl ) equal to 0.9 and 0.5 respectively; (d) error treatment has been applied to all 
detections (A(lg 10 iVjji) = 0). Case (c) has been used for the results derived in Paper I. Fig. 3 
shows the resulting value of a for the 4 cases and the comparison between the theoretical 
and observational cumulatives. Notice that, when the treatment of errors is operative, there 
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is a good match between the two cumulatives. The "observational" cumulatives in these 
cases (dotted histograms, in Fig. 3) depend on the assumed distribution, because the A^jji 
positions vary accordingly. 

The logarithm of the Likelihood C a , normalized to its maximum value, is shown in 
Table 1 as we vary a for all the 4 cases. We can see that the maximum of the Likelihood is 
very well determined, without the presence of secondary maxima, and with 1-a uncertainties 
of the order of 0.1. In Table 1 we also give the significances Cr and Cc of the uniformity 
of the Rj and Cj distributions. Rj and Cj can be uniformly distributed for a very wide 
range of a when the number of detections with variable iVjjj is large; however only near the 
Maximum Likelihood value of a Rj and Cj are both uniformly distributed. The maximum 
of Cc is better defined than that of Cr, and its position agrees more closely with that of the 
maximum value of C Q . The uniformity of the Cj distribution represents a valid test of the 
goodness of the Maximum Likelihood solution and of the fitting routine. 

In Paper I we have shown that, in the a-X plane, the constant-probability contours are 
elongated, giving a global uncertainty on a of about ±0.3 when X uncertainties are taken 
into account. Fig. 4 shows the probability contours in the a-X plane for 3 redshift intervals 
which together cover all the redshifts in our compilation, each containing the same number 
of detections (~ 90). We find that at lower redshifts (Subsample A) the probability contours 
are shifted towards higher values of a and X; at higher redshifts (Subsample C) instead the 
allowed a values do not show significant changes, while X decreases. However, since these 
variations are of the same order of parameter uncertainties, they could just be statistical 
fluctuations. Further insight will come out from the simulated data analysis in the next 
session. 

An issue that will be discussed more in detail in a future paper is whether line blending 
is present in the data and could affect the above results. The probability of blending of 
Lyman-a lines is relatively low for randomly distributed absorbers, but clustering or blending 
with some metal lines might imply a somewhat higher probability and affect lines with 5 < 
W < 10 A, sometimes declared as non damped. To address this point we have performed an 
analysis of the data without considering low-resolution searches for lines with 5 < W < 10 A: 
the a and X values that we find are consistent with those relative to the whole data sample. 
Therefore we tend to believe that our data sample closely reflects the real Hi distribution. 
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18 19 20 21 18 19 20 21 

Fig. 3. — The Maximum Likelihood values of a and the corresponding cumulative number 
of absorbers in the z range [1.75,3.25] are shown for different prescriptions of the errors 
treatment. The solid line is the theoretical cumulative function, while the dotted line is the 
observational one. 




Fig. 4. — 1-cr, 2-cr and 3-cr confidence levels in the X-a plane, for 3 z-selected subsamples, 
labeled by A (dotted lines), B (solid lines) and C (dot-dashed lines). 
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5. Simulations 

Possible biases might affect the data analysis, either related to our assumptions about 
the physical properties of the absorbers, or to our interpretation of the observations. In this 
section we show how to use simulated data for investigating the latter type of biases. These 
might derive from: (i) an uneven coverage of iVjjj and z since the database contains results 
from surveys with different sensitivities; (ii) non uniform uncertainties in iVjjj; (in) incorrect 
estimates of detection thresholds or of uncertainties. Simulated data can also be used 
to test the statistical procedure described in the previous section, and the confidence levels 
properties for the best fitting parameters. 

5.1. The numerical procedure 

Simulations supply, in a controlled way, a large number of samples to which we can 
apply our statistical analysis. We have devised a routine to simulate data in a highly real- 
istic way: starting from an "a priori" distribution function for the absorbers total column 
densities, it generates a data set compatible with a given column density and redshift cov- 
erage. In the present work we use a coverage consistent with the observations listed in our 
database. The basic ingredients of the simulation are: 1. an original distribution function 
of -/Vji_|_ for absorbers of a given geometrical shape; 2. an iVjjj_-iVjjij_ relationship for each 
redshift bin; 3. a N^-z coverage for the simulated survey. After the total face-on column 
density has been transformed into a line of sight Hi value, we introduce measurement errors 
similar to those present in the actual survey compilation. These errors are simulated by 
deriving first the random quantities Rj, uniformly distributed, and then locating the uncer- 
tainties [Nn ljjnin j, iVjji moa ; j] such that the position of "true" values iVfjj corresponds to Rj. 
When (lg 10 N-gj min j + lg 10 Ny l max j)/2, is above 10 17 - 68 cm -2 the absorption is considered as 
saturated. In this case, if there has been a search for the corresponding damped Lyman-oj 
but the line has not been found, the uncertainty range is set to [10 17 ' 68 , 10 20 ' 13 ] cm" 2 ; if in- 
stead the damped Lyman-a line has not been searched for, the uncertainty range is set to 
[10 17 - 68 , 10 21 - 7 ] cm" 2 . 

For saturated absorption we simulate also the shadowing effect. Shadowed absorbers are 
eliminated from the sample and the N-gj-z coverage of the survey is modified accordingly for 
LLS and DLS searches: this "after-shadowing" coverage is what should be compared with 
the real data coverage. 
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5.2. Simulation results 

We first simulate samples of data consistent with the Maximum Likelihood solution 
a = 2.683, X = 2.82, obtained for real data in the z range [1.75,3.25] for pi = 0. Results 
of the statistical analysis on one of our simulated data samples are shown in Fig. 5, where 
the efficiency of our method in recovering the original distribution is clear. Fig. 5 is the 
analogous of Fig. 3 for real data, the only difference being that now in (a) the "true" values 
of iVjji for simulated data are used, without any degradation by measurement processes and 
are compared directly with the model. The slight disagreement between the two curves in 
(a) is indicative of the magnitude of statistical fluctuations in the simulation. 

C , Cr and Cc for the simulated sample show a behavior very similar to that for real 
data (Table 1). An increase in the number of detections on which the treatment of errors 
is applied improves the match between the "measured" distributions and the original ones. 
Our algorithm is very efficient in recovering the original Rj distribution as well but there is 
no way to recover the original Rj values of individual detections. Detections with the same 
range of column densities are in all respects identical among themselves, and any permutation 
between their respective Rj values cannot be traced (i.e. there is no correlation between the 
"true" and derived Rj values). 

To investigate if the selection of an error threshold affects the determination of a we 
have produced 5000 different simulated samples with a = 2.683. On each sample we apply 
the analysis with the 4 different error thresholds as in Fig. 5. The results, reported in 
Table 2, show that the best fit a values are only slightly smaller than the original ones, and 
the dispersion for the best a values found (cr(cx)) is consistent with the average estimated 
uncertainty (< cr >). 

We can also use the simulated data to check the uncertainties in the best-fitting pa- 
rameters (cx,X) by producing a confidence level map for the simulated data, for 3 ^-selected 
subsamples, analogously to what has been done in Section 4.3 for the real data (Fig. 4). 
Fig. 6 shows for one simulated sample that the confidence levels present an elongation simi- 
lar to the levels in Fig. 4. This means that in the real data there is no evidence of deviations 
of the gas distribution from our model assumptions. The elongation depends on a partial 
degeneracy of models in the a-X plane: a should increase as X increases to give the same 
slope for the distribution of low Hi column density LLS, a drawback that will be difficult to 
cure by a moderate increase of the numbers of known absorbers. 

By looking at the three subsamples (Fig. 6) it is clear that statistical fluctuations are 
present mostly in the direction of the elongation, consistently with the estimated uncertain- 
ties. With the present sample coverage the real data are still consistent with no evolution of 
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(a) J: 




I I I I I I I I I I I I I I I I I I I I I I 



(c) J: 




18 19 20 21 18 19 20 21 



Fig. 5. — Various fits to the cumulative number of absorbers in the z range [1.75,3.25], on 
simulated data. This figure is analogous to Fig. 3, except for panel (a), which gives the 
"true" cumulative distributions, both theoretical (solid line) and "measured" (dotted line). 




Fig. 6. — Map of the confidence level in the X-a plane, relative to a simulated sample, for 
3 ^-selected subsamples. All definitions in this figure are analogous to those in Fig. 4. 
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a and X, but some evolution cannot be excluded, especially for the high redshift sample. For 
this sample in fact the shift of the probability contours is not in the direction of elongation 
but mostly towards lower values of X. This would imply that either the ionizing flux is lower 
or that the absorbing gas has higher volume densities than at lower redshifts. Our analysis 
shows that in order to have a definite answer on the evolution of a and X closer contour 
levels are needed. This can be achieved by determining observationally the column density 
of saturated LLS (through lines in other spectral regions or through the recover of the quasar 
flux at shorter wavelengths) or by increasing the number of known high iVjjj DLS. 



-23- 



REFERENCES 

Auriere, M. 1982, A&A, 109, 301 

Bahcall, J. N., Bergeron, J., Boksenberg, A., Hartig, G. F., Jannuzi, B. T., Kirhakos, S., 
Sargent, W. L. W., Savage, B. D., Schneider, D. P., Turnshek, D. A., Weymann, R. 
J., and Wolfe, A. M. 1993, ApJS, 87, 1 

Bechtold, J., Green, R. F., Weymann, R. J., Schmidt, M., Estabrook, F. B., Scherman, R. 
D., Wahlquist, H. D., and Heckman, T. M. 1984, ApJ, 281, 76 

Canizares, C. R., Grindlay, J. E., Hiltner, W. A., Liller, W., and McClintock, J. E. 1978, 
ApJ, 224, 39 

Corbelli, E., Salpeter, E.E., and Bandiera, R. 2001, ApJ, 549, in press (Paper I) 

Courvoisier, T. J.-L., and Paltani, S., eds. 1992, IUE-ULDA Access Guide No. 4A (ESA SP 
1153) 

Cristiani, S., D'Odorico, S., D'Odorico, V., Fontana, A., Giallongo, E., and Savaglio, S. 1997, 
MNRAS, 285, 209 

Djorgovski, S., and King, I. R. 1984, ApJ, 277, L49 

Hagiwara, K., and Zeppenfeld, D. 1986, Nucl.Phys., 274, 1 

Harris, W. E., and van den Bergh, S. 1984, A J, 89, 1816 

Henon, M. 1961, Ann.d'Ap., 24, 369 

Jannuzi B. T., Bahcall, J. N., Bergeron, J., Boksenberg, A., Hartig, G. F., Kirhakos, S., 
Sargent, W. L. W., Savage, B. D., Schneider, D. P., Turnshek, D. A., Weymann, R. 
J., and Wolfe, A. M. 1998, ApJS, 118, 1 

Kauffmann G. 1996, MNRAS, 281, 475 

Kendall, M.G., and Stuart, A. 1969, "The Advanced Theory of Statistics" Vol. 2 (Griffin: 
London) 

King, I. R. 1966, AJ, 71, 276 

King, I. R. 1975, Dynamics of Stellar Systems, A. Hayli, Dordrecht: Reidel, 1975, 99 
King, I. R., Hedemann, E., Hodge, S. M., and White, R. E. 1968, AJ, 73, 456 



-24- 



Kron, G. E., Hewitt, A. V., and Wasserman, L. H. 1984, PASP, 96, 198 
Lanzetta, K. M. 1988, ApJ, 332, 96 
Lanzetta, K. M. 1991, ApJ, 375, 1 

Lanzetta, K. M., Wolfe, A. M., Turnshek, D. A., Lu, L., McMahon, R. G., and Hazard, C. 
1991, ApJS, 77, 1 

Lanzetta, K. M., Wolfe, A. M., and Turnshek, D. A. 1995, ApJ, 440, 435 

Lynden-Bell, D., and Wood, R. 1968, MNRAS, 138, 495 

Newell, E. B., and O'Neil, E. J. 1978, ApJS, 37, 27 

Ortolani, S., Rosino, L., and Sandage, A. 1985, AJ, 90, 473 

Peterson, C. J. 1976, AJ, 81, 617 

Rauch, M. 1998, ARA&A, 36, 267 

Sargent, W. L. W., Steidel, C. C, and Boksenberg, A. 1989, ApJS, 69, 703 

Savaglio, S., Cristiani, S., D'Odorico, S., Fontana, A., Giallongo, E., and Molaro, P. 1997, 
A&A, 318, 347 

Spitzer, L. 1985, Dynamics of Star Clusters, J. Goodman and P. Hut, Dordrecht: Reidel, 
109 

Stengler-Larrea, E. A., Boksenberg, A., Steidel, C. C, Sargent, W. L. W., Bahcall, J. N., 
Bergeron, J., Hartig, G. F., Jannuzi, B. T., Kirhakos, S., Savage, B. D., Schneider, 
D. P., Turnshek, D. A., and Weymann, R. J. 1995, ApJ, 444, 64 

Storrie-Lombardi, L. J., McMahon, R. G., Irwin, M. J., and Hazard, C. 1994, ApJ, 427, L13 

Storrie-Lombardi, L. J., McMahon, R. G., Irwin, M. J., and Hazard, C. 1996, ApJ, 468, 121 

Storrie-Lombardi, L. J., Irwin, M. J., and McMahon, R. G. 1996, MNRAS, 282, 1330 

Storrie-Lombardi, L. J., and Wolfe, A. M. 2000, ApJ, 543, 552 

Turnshek, D. A., Wolfe, A. M., Lanzetta, K. M., Briggs, F. H., Cohen, R. D., Foltz, C. B., 
Smith, H. E., and Wilkes, B. J. 1989, ApJ, 344, 567 

Tytler, D. 1982, Nature, 298, 427 



-25 - 

Tytler, D. 1987, ApJ, 321, 49 

Weinberg, D. H., Miralda-Escude, J., Hernquist, L., and Katz, N. 1997, ApJ, 490, 564 
Wolfe, A. M., Turnshek, D. A., Smith, H. E., and Cohen, R. D. 1986, ApJS, 61, 249 
Wolfe, A. M., Lanzetta, K. M., Foltz, C. B., and Chaffee, F. H. 1995, ApJ, 454, 698 



This preprint was prepared with the AAS IATgX macros v5.0. 



-26- 



A. Averaged cross section for homogeneous ellipsoids 

In this Appendix we compute analytically the differential cross section for absorbers 
which can be modelled as axisymmetric ellipsoids of constant density. We use a cartesian 
coordinate system (x,y,z) in which z is parallel to the line of sight. The ellipsoid has a 
semi-minor axis h and semi-major axes R along the orthogonal directions. If its axis of 
symmetry is tilted along the y axis by an angle 6 with respect to the line of sight, the surface 
of the ellipsoid is: 

/sin 2 # cos 2 #\ 2 n . „ n (\ 1\ /cos 2 # sin 2 #\ , 1 y 2 , Ail , 

+ ^)x 2 + 2^0cos9 -- — )xz+ — + — )z 2 = l-^-. (Al) 



V h 2 R 2 J \h 2 R 2 J \ h 2 R 2 J R 2 

The spatial extent of the ellipsoid in the z direction is 2(. Defining: 

£ 2 = cos 2 + sin 2 0-^, (A2) 



e=tU-^L-t)- (A3) 



we can write: 

t 2 = 

i 2 v i 2 R 2 R 2 , 

The column density along the line of sight is A^jjj = 2n(. Since ( < h/C, the distribution 
for a fixed orientation presents an upper cut of A^jjij_/C (where Nn l± = 2nh is the maximum 
face-on column density). 

The ellipsoid cross section for column densities larger than Njjj is: 

e( ^^H^) 2 ) <A4) 

which gives the differential cross section for one fixed orientation in space: 

a(N m ,») = nR 2 ^, (A5) 
(yU = cos 6) and an orientation-averaged cross section: 

a(N m ) = 7i R 2 f 1 6 (N m± - £7V Hl ) dp. (A6) 



'0 



Introducing the variables: 



« = T = t + {w- 1 )* H/ =f mta ( 1 'lf | - (A7) 
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the average differential cross section for randomly oriented absorbers is: 

, x _ 2h 3 N m f w w 4 

a(N m ) = 71 R =^^= / dw 

1 mj N^ l± R 3 ^R 2 /h 2 -lJi 

3h A N m [W ( 1 + -W 2 J VW 2 - 1 + lg (V + V^ 2 - l) J 
= ttR 2 ^ 1 (A8) 

For Ajjj < 2n/i the ^(iVjjj) increases linearly with A" Hl while it vanishes for > 2nR. In 
the intermediate range the behavior of <x is more complex, but for Ajjj <C 2nR, the following 
approximation holds: 

HNm) = I 1 - -CM | --Hi- j II. (A9) 




/V 3 

where A = 7ri? 2 /2y / l — h 2 /R 2 . Notice the similarity between Eq. (A9) and Eq. (6) for 
infinitely thin slabs. In Fig. 7 we compare the exact solution for the ellipsoids of various 
axial ratios h/R to the approximate formula we use in our routines (Eq. (7)). 

Although a depends on the detailed shape of the absorbers, the approximation by a 
power law with an upper cutoff is appropriate for ellipsoids as well as for other geometrical 
shapes, especially for thickness h/R <^ 1. The average value of the projected column density 
is instead more geometry dependent. For ellipsoidal absorbers this is: 



Io m±R/h N m a(N m )dN m 
J N ^ R/h a(N m )dN m 

4N m± RVR 2 - h 2 



(AlO) 



3 [r^/R 2 -h 2 + h 2 \n^R + VR 2 - h 2 } /h^y 
For h <C R, Nui approaches 4Ajj i± /3. This result differs from that for a slab, for which: 

A^Hi = 2AW (1 " h/R) j (1 - h 2 /R 2 ) , (All) 
whose limiting value is N^i = 2A / h i _ l for infinitesimal thickness. 
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lgio(N„,/N HIf ) 



Fig. 7. — Comparison between orientation-averaged cross sections for ellipsoids (solid curves) 
and those for the finite slab approximation (dotted curves), for various axial ratios (h/R = 
0.3,0.1,0.03,0.01). 
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Table 1. Best fits to the data a 





Odoc la J 




Case (b) 






Case (c) 






Case (d) 






A II n/^l T1TC flVO/l 

puiiiis n.x.eu. 




?10 7V HI ) > .9 


A He 


5io^m)>0.5 


All points var. 


cx 


In r 

in L, 


In r 

111 L, 


In C R 


\nC c 


In r 




ln£ c 


ln£ G 


In C R 


\nC c 


2.05 


-25.2 


-21.5 


-0.083 


-6.218 


-19.7 


-0.022 


-0.784 


-18.1 


0.000 


0.000 


2.15 


-18.1 


-15.0 


-0.078 


-3.784 


-13.7 


-0.014 


-0.253 


-12.4 


0.000 


0.000 


2.25 


-12.3 


-9.91 


-0.066 


-1.999 


-8.90 


-0.003 


-0.039 


-7.94 


0.000 


0.000 


2.35 


-7.73 


-5.94 


-0.034 


-0.939 


-5.25 


0.000 


-0.001 


-4.55 


0.000 


0.000 


2.45 


-4.31 


-3.05 


-0.011 


-0.311 


-2.61 


0.000 


0.000 


-2.15 


0.000 


0.000 


2.55 


-1.93 


-1.16 


-0.002 


-0.054 


-0.92 


0.000 


0.000 


-0.67 


0.000 


0.000 


2.65 


-0.52 


-0.18 


0.000 


-0.006 


-0.11 


0.000 


0.000 


-0.03 


0.000 


0.000 


2.75 


0.00 


-0.05 


0.000 


-0.003 


-0.09 


0.000 


0.000 


-0.18 


0.000 


0.000 


2.85 


-0.31 


-0.71 


0.000 


-0.006 


-0.82 


0.000 


0.000 


-1.04 


0.000 


0.000 


2.95 


-1.39 


-2.10 


0.000 


-0.025 


-2.24 


0.000 


0.000 


-2.57 


0.000 


0.000 


3.05 


-3.20 


-4.18 


0.000 


-0.225 


-4.31 


0.000 


0.000 


-4.73 


0.000 


0.000 


3.15 


-5.68 


-6.90 


0.000 


-0.851 


-6.98 


0.000 


-0.015 


-7.47 


0.000 


0.000 


3.25 


-8.80 


-10.2 


0.000 


-1.857 


-10.2 


-0.002 


-0.138 


-10.8 


0.000 


-0.008 


3.35 


-12.5 


-14.1 


-0.000 


-3.175 


-14.0 


-0.005 


-0.476 


-14.6 


-0.001 


-0.131 


3.45 


-16.8 


-18.5 


-0.001 


-4.780 


-18.3 


-0.014 


-1.067 


-18.8 


-0.005 


-0.504 


3.55 


-21.6 


-23.4 


-0.001 


-6.648 


-23.0 


-0.032 


-1.929 


-23.6 


-0.010 


-1.122 



Fits for X = 2.82, h/R = 0.2, 7 = 1.0, lg 10 N th = 17.20 and data in the redshift range [1.75,3.25]. 
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Table 2. Results for 5000 simulations 



Case 


a 




a(a) 


< a > 


(a) 


2.549 


-0.134 


0.116 


0.106 


(b) 


2.468 


-0.215 


0.113 


0.106 


(c) 


2.569 


-0.114 


0.122 


0.116 


(d) 


2.554 


-0.129 


0.120 


0.116 



